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Abstract. The problem of resolving spherical harmonic components from numerical 
data defined on a rectangular grid has many applications, particularly for the problem 
of gravitational radiation extraction. A novel method due to Misner improves on 
traditional techniques by avoiding the need to cover the sphere with a coordinate 
system appropriate to the grid geometry. This paper will discuss Misner's method and 
suggest how it can be improved by exploiting local regression techniques. 
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1. Introduction 

Spherical harmonic decomposition is an important tool in computational science. In 
numerical astrophysics, for example, harmonic decomposition can serve as the basis 
of an algorithm for extracting gravitational waveforms from simulations of radiating 



Harmonic decomposition of data defined on a two-sphere is a straightforward 
computation. Consider a scalar field $(t, r , 0, 0) defined on a two-sphere of radius 
r . The field can be expanded in terms of scalar spherical harmonics as 



£,m 

The orthonormality of the spherical harmonics with respect to integration over the 
sphere can be used to compute the harmonic amplitudes 



where the bar denotes complex conjugation. 

Of interest for numerical analysis are cases in which the field <3> is not known as 
an analytic function on the sphere, but rather as the solution of, for example, a wave 
equation evolved on a discrete set of points {xi} which comprise a three-dimensional 
grid for the numerical simulation. In general these grid points will not lie on the sphere 
r = ro. Therefore, computing the harmonic coefficients <E> £m from Q requires a method 
of determining field values on the integration sphere. 

The most obvious method of computing <3> £m from grid data is first to interpolate 
data from the grid onto coordinate patches covering the two-sphere, then compute a 
numerical integral of the interpolated data against the conjugate spherical harmonics. 
Because the grid data is the source of all information about $, it is important that the 
coordinate patches on the two-sphere adequately represent the distribution of three- 
dimensional grid points near the sphere. For example, a spherical coordinate grid 
discretized in {8, <ft} will over-represent points near the poles and under-represent points 
near the equator. 

Previous investigations of the the gravitational radiation extraction problem used 
a pair of stereographic coordinate patches to cover the sphere p. These patches, due to 
Gomez, et al jl], overlapped at the equator and special care was required to compute 
the integrals of ((H) correctly. A useful alternative is to use "cubed sphere" coordinates 
like those used by Zink, et al In these coordinates the six patches represent 

cartesian grids, deformed to cover a sphere, which join nicely at the patch boundaries. 
The distribution of points in these patches is approximately uniform and close to the 
distribution of nearby grid points. 

Misner has proposed an attractive alternative approach to the multipole 
decomposition problem that replaces the traditional two-step interpolation/integration 
method with a single volume integration step for each mode [3]. In addition to the 
computational simplification offered by this approach, Misner's method completely 



systems. HJEIE! 



$(t,r„,M) = ^^ m (t,r )^ m (M) • 




(1) 
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bypasses the need to choose a coordinate grid covering the two-sphere. Instead, it 
relies directly on values of $ on the grid points surrounding the sphere. The relative 
weight of each point in the calculation depends only on its proximity to the sphere. 

Fiske jH] has recently presented an analysis of symmetry and convergence properties 
of Misner's method, and Fiske, et al have demonstrated its usefulness in extracting 
gravitational radiation from cubical grids with fixed mesh refinement. 

This paper presents an alternative analysis of Misner's algorithm which exploits 
the fact that it can be cast as a weighted local least-squares regression procedure. 
This approach somewhat simplifies the computational aspects of Misner's method, and 
improves the accuracy of the calculations. 

While the analysis in this paper will assume complex-valued harmonics, test 
computations on a real-valued scalar field will be performed using real-valued spherical 
harmonics for simplicity. The test function will be the analytical quadrupole solution 
of a three-dimensional linear scalar wave equation. Specifically, the test function will 
be <&(t,r, 8, <fr) = Q 20 (t,r)Y2o(6,<j)), where $ 20 is constructed from a generating function 
of the form (x/X 2 ) exp(— x 2 /X 2 ), and x = t ± r. The initial data is time-symmetric and 
we choose unit wavelength (A = 1) and amplitude. This is similar to the wave evolved 
by Fiske, et al j2j, but here $ is the solution of a linear, scalar wave equation. The 
computational grid extends from —8 to 8 with uniform spacing h in each dimension. 
Multipole amplitudes are computed from this test data on a two-sphere of radius ro = 6. 
Most computations are performed in the first octant with even reflection symmetry 
imposed at the coordinate plane boundaries. 

2. Misner's Algorithm 

2.1. Continuum Limit 

Recognizing that a surface integral like can be treated as the derivative of a 
volume integral, Misner approaches the problem of computing harmonic coefficients 
by integrating data over a specified volume surrounding the two-sphere. The volume of 
interest is a spherical shell, S, of half-width A that surrounds the sphere: 



He then expands the field $ not only in terms of orthonormal basis functions 
in {#,</>} (spherical harmonics), but also in terms of radial basis functions R n {r) 
(n = 0, . . . , oo) that are orthonormal with respect to integration over the radial interval 



The tensor product of the radial and angular basis functions represents a new set of 
three-dimensional basis functions that are orthonormal with respect to integration over 
the shell S. These can be written as 



S = {x = {r,9,4>} \ r e [r - A,r + A]} . 



[r - A,r + A]: 




(2) 



Y A = Ynl 



m 



Rn{r)Y tm {e,(t>) , 
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where the single index A = {n£m} denotes the three indices that identify the three- 
dimensional basis functions. 

Specifically, Misner uses Legendre polynomials P n to generate the (real-valued) 
radial basis functions 



r V 2A " V A 

which satisfy (j2J). The choice of Legendre polynomial as the basis for R n is not unique 
and other choices will be investigated later. 

The integral over S defines an inner product for functions defined on the shell: 

(f\g) = [ f(x)g(x)d 3 x, 
Js 

where the bar denotes complex conjugation. The ort honor mality relation can be 
expressed as (Ya>\Y a ) = 5 A > a = 5 n > n 8t>i8 m > m . 

The expansion of $ in terms of these basis functions is 

$(t,x) = ^$ A (t)r A (r,M), 

A 

where the coefficients of the expansion are <& A (t) = (Ya|$). The harmonic coefficients 
at r are expressed in terms of the coefficients $ A as: 

$ £m (t,r ) = J2 R n(r )® nem (t) = Y, R n(ro) (Y ntm \®) • (3) 

n n 

It is useful to re-cast (j3J) in terms of a set of projection functions, 



'tm(x) , (4) 

that act on $ under the inner product to give the harmonic amplitudes: 

$ £m (t,r ) = (j/ m |$) = I p em (x)<f>(t,x)d 3 x . (5) 

Js 

2.2. Numerical Approximation 

Misner's algorithm is based on a specific method of converting the continuum integral 
above into a numerical integral over the grid points located within the shell S. By 
assuming a cell-centered grid of uniform spacing h (cell volume h 3 ), one can approximate 
a volume integral by a sum that is equivalent to a "midpoint" quadrature: 

Rxi) g(xi) /i 3 w f(x)g(x)d 3 x, (6) 

ignoring effects near the boundary of S. 

Of course some points near the shell boundary belong to cells that do not lie entirely 
within S, while other points located just outside the boundary belong to cells that lie 
partially within the shell. To account for this effect, Misner adds to the sum © all 
points a distance less than \h outside of S. The correct inner product uses points from 
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a larger shell, which we shall denote S + , defined to be the set of all grid points that lie 
within a radial distance 5 = A + |/i of ro. Furthermore, all points that lie "near" the 
shell boundary are weighted based on the partial volume of their cell that lies within S. 

Keeping track of these partial volumes is complicated, but a simple approximation 
is to treat each point near the boundary as if it lay on a coordinate axis. This makes 
the partial volume a linear function of r: 

{0 \n - r | > 5 

h 3 \vi — r | < 5 — h 

(5 — \ri — ro|) h 2 otherwise, 

where is the radial coordinate of point Xj. Other choices of weight function will be 
discussed below. 

In Misner's method, computations of the harmonic coefficients from grid data are 
based on the numerical inner product 

(f\g) = f(xi)g(xi)w(n). (7) 

In the continuum limit, the inner product of basis functions is equivalent to the 
identity matrix. However, because the weighted sum of (J7J) is only an approximation 
to the continuum integral, the basis functions Ya are not strictly orthonormal and the 
numerical inner product defines a matrix Gab = (Ya\Yb) that approximates the identity. 
While some matrix elements of Gab will be identically zero due to symmetries, the 
remaining elements will only approximate their continuum values to a degree determined 
by the numerical approximation. As a result, using the numerical inner product (J7J) in 
the calculation of (jHJ) gives a convergent set of values of $^ m (t, r ), but could result in 
mode mixing. 

We can avoid this by using the inverse of the inner product matrix, G AB = {Gab) 1 , 
to define a set of dual basis functions Y A = G AB Yb that are orthonormal to the 
original basis with respect to the inner product: 

(Y A \Y B ) = J2G AC (Y C \Y B ) = Y,G AC G CB = S A . 

c c 

Using the dual basis functions, define the expansion coefficients $ using the inner 
product with the dual basis 

$-4 = (y a \$) = J2 gAB ( Y b\®) , (8) 

B 

using the numerical inner product This choice defines an approximation to the 
original data, $ = ^ A <& a Ya ~ 3>. For a given finite collection of basis functions, 
the approximation $ minimizes the (weighted) sum of the squared differences between 
and over the grid points in S + . In other words, Misner has replaced the 

interpolation/spherical integration with a least squares procedure to compute spherical 
harmonic coefficients. 

The projection functions of (@J) can be computed at the beginning of a simulation 
from the dual basis functions Y A and stored for future use. At each time step of a 
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simulation, the approximate spherical harmonic amplitudes $^ m (t,r ) are computed 
from the numerical inner product of the projection functions with the field <3? 

$ £m (t,r ) = (j/ m |$) = ^(xjQ&xJwin) . 

Xi£S + 

3. The Least-squares Perspective 

In their analyses, Misner and Fiske focus moslty on the relation between the weighted 
sums and the underlying three-dimensional integrals they approximate. It is useful, 
however, to consider Misner's method strictly as a least-squares problem. 

3.1. Matrix Form 

Focusing on the least-squares aspect of the method makes it easier to see how to 
compute the projection operators that will be used in the computation of the harmonic 
amplitudes. It is useful to re-cast the algorithm in the traditional matrix notation of 
linear regression analysis. 

Let N be the number of points in the computational shell S + , and M be the 
number of independent modes A = {n, £, m} over which the three-dimensional harmonic 
expansion is performed. Let ^ be the N x 1 column vector of values of <3> at all points 
in S + . Similarly, let Y be the N x M design matrix, where each column of Y is one 
of the M basis functions Ya evaluated at all points in the shell. Finally, let W be the 
N x N diagonal matrix of weights u>.j. Then Gab is easily computed as G = Y'WY, 
where Y* is the Hermitian conjugate of Y . 

Denote the set of expansion coefficients $ A as the (M x 1) column vector F. 
The vector of approximate values, $ ~ $ in S + , is & = YF. The weighted sum 
of squared residuals (SSR) between $ and $ is (3> — YF)' W (<P — YF), and the 
expansion coefficient vector that minimizes the SSR is the familiar matrix solution of 
the least squares problem, 

F = (Y*WY)~ 1 Y*W& = G l Y ] W$ . (9) 

This is obviously equivalent to (JBJ) above. 

Let P = G 1 Y^W be the projection matrix that computes the expansion 
coefficients when multiplied by the grid data: F = P As in (J3J), appropriate 
linear combinations of rows of P define a set of row vectors p lm that, when multiplied 
by the grid data, give the spherical harmonic amplitudes $^ m (t, r ) = p £m <I>. These 
projection vectors are computed at the beginning of a simulation and stored for later 
use in computing the time-dependent expansion coefficients as the simulation evolves. 

Of course it is not generally advisable to compute P by inverting G. It is better 
to compute the M x N matrix P as the solution of the linear system 

GP = Y ] W . 

Because G is a real, symmetric, positive-definite matrix, this system can be solved 
rapidly using Cholesky decomposition. 
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S 1 .^. Weighting and Local Regression 

There are actually two types of least-squares regressions involved in the method: a 
global regression of all points over the spherical harmonic basis functions, and a local 
regression over the radial basis functions. Local regressions differ from traditional global 
regression fits in that they are only used to compute the best-fit function at a specific 
point. In this case the local regression computes the spherical harmonic amplitudes at 
r = To from a fit over the interval (ro — 5, ro + 5). 

Local regression techniques are popular for data smoothing algorithms [7] . In such 
cases the value of a data point is replaced by the result of a weighted polynomial fit to 
it and neighboring points. A common, though not universal, feature of local regression 
techniques is that the weight of each point in the fit decreases with distance from 
the evaluation point. Two popular choices of local weight function are the "bisquare" 
function, 

B(u) = 

and the "tricube" function 
T(u) = 

where u = (r — ro)/5 is a dimensionless radial coordinate defined so that the domain 
(r — 5, r + 5) is transformed to the interval (—1,1). 

Contrast this with Misner's method. In the continuum limit, all points weigh 
equally in the computation of G because the Legendre polynomials are orthonormal 
with respect to integration over a constant weight function. 

Now consider Misner's weight in the numerical approximation, written in terms of 
the dimensionless radial coordinate u. Up to a multiplicative factor of h 3 the weight 
function is 

{0 \u\ > 1 

1 \u\ < 1 - h/5 

(1 — \u\) 5/h otherwise. 

While the weight is constant for points near r , it falls off linearly for points near the 
edges of the shell S+. The degree of down- weighting depends on the ratio of the size of 
the shell to the grid spacing, 5/h. For 5 ^> h, the weights are essentially constant. 

Figure Q shows Misner's weight function for two different values of 5/h and 
compared to the bisquare and tricube weights. In his analysis of Misner's method [H], 
Fiske finds that 5 = jh (equivalently, A = jh) gives good results in tests. In this 
case most points in S+ are "near" the boundary, and the weight function is local in the 
sense that it only emphasizes points nearest the center. In fact, Fiske's choice leads to a 
weight function that is very similar in profile to the traditional local regression weights, 
and gives similar results. 
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Figure 1. Comparison of Misncr's weights and local weights. The functions are scaled 
to have the same area underneath each curve. 



To see the benefit of using local weights, consider Fiske's analysis in the continuum 
limit. Misner's method works because a linear combination of the radial basis functions 
R n form a Dirac delta function: 

oo 

J2 Rn(r ) Rn(r) = r- 2 5(r - r) . 

n=0 

Given the definition of R n in terms of the Legendre polynomials, define 

n=0 ^ ' 

Note that only even values of n figure into this sum because P n (0) = for odd n. The 
function dp{u,N) will approach the Dirac delta function S(u) in the limit as iV — > oo. 
The use of a finite number of radial basis functions produces a truncation error that will 
have leading term 0(5 N+2 ), so that a fit with iV = 2 or 3 will converge with shell size 
as 0(5 4 ). 

The same analysis can be applied using functions that are orthonormal with 
respect to integration over the bisquare or tricube weight functions. Such orthogonal 
polynomials can be constructed using a Gram-Schmidt procedure. At least in the case 
of bisquare-weight functions, they can also be found by the simple Rodrigues formula 

M«) = c B (i-« 2 )- 2 ^(i-« 2 ) B+2 , 

where c„ is a normalization constant. 
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Figure 2. Comparison of d(u, N) constructed from Legendre and bisquare-weight 
polynomials. 



The first few normalized bisquare-weight polynomials are 

b (u) = J% h(u) = J^u 



b 2 (u) = y/% (1 - 7m 2 ) b 3 (u) = (u - 3u 

h{u) = (1 - 18m 2 + 33m 4 ) . 



Define 



N 



d b {u,N) = J2b n {0) K 



[u] (1 — u 2 ) 2 



n=0 



Again, this will approach S(u) in the limit N ^ oo and the leading order of truncation 
error is still 0(5 N+2 ). 

Figure 121 compares dp and df, for low orders. It is clear from the figure that the 
bisquare approximation to 8{u) is better at each order than the Legendre approximation. 
This is because of the down-weighting inherent in local regressions. 

The advantage of the local weighting is apparent in numerical tests. Figure El 
compares results of computations of $ 20 from the test function using Misner's weight 
to those using the tricube weight. (Errors using the bisquare weight are close to, but 
slightly larger than, those using the tricube weight.) The vertical axis is the (base- 10) 
log of the L 2 norm of the error, while the horizontal axis is the log of S/h. The tests 
were performed using a quadratic fit in r (N = 2). 

The error in both plots in figure El is fourth-order convergent with respect to 5 
(for fixed h), as predicted by Fiske's analysis. The error using Misner's weight actually 
converges slightly faster than fourth-order because the shape of Misner's weight function 
changes with S/h, while the shape of the tricube weight is fixed. As S decreases, Misner's 



Harmonic Ampitudes From Grid Data 



10 



-1 



Misner's Weight 

Tricube Weight 




-2 ■ 



-3 ■ 



-4 ■ 



-0.2 







0.2 



0.4 



0.6 



Figure 3. Comparison of log(error) to \og(S/h) for h = 1/8. The tricube error is 
4th-order convergent, and the error for Misner's weight converges slightly faster. 



weight function becomes more local, and it should be no surprise that Fiske found good 
results with 8 = Regardless of the value of S/h, however, the bisquare and tricube 
weights will give better results than Misner's weight. 

3.3. Switching to the monomial basis 

The preceding analysis using Legendre and bisquare-weight polynomials is strictly 
correct in the continuum limit, not for calculations from grid data. In fact, the properties 
of Legendre polynomials in the continuum limit cannot be expected to carry over into 
the numerical regime, especially for small values of S, because the weight function will 
no longer be constant. Furthermore, because Misner's down-weighting prescription for 
boundary points is only correct for points near the coordinate axes, the orthonormality 
properties of the Legendre polynomials will break down in the numerical approximation 
for small shell sizes. 

Of course this is not an problem in practice because Misner's method doesn't really 
rely on the computation of integrals. Instead, it is based on solving a weighted least 
squares problem and the only requirement for the radial basis functions is that they 
be linearly independent. It is quite simple (and useful in practice) to switch from an 
orthogonal polynomial basis to a monomial basis: 



This choice simplifies the calculation by removing the need to form a linear combination 
of rows of P to get the projection vectors p em . On the contrary, because Ro(r ) = 1 
and R n (r ) — for n > 1, then $ £m = $ 0£m and the rows of P corresponding to n = 




(10) 
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Figure 4. Comparison of errors in computation of d$ 20 /dr for quadratic and cubic 
fits. The plot shows log(error) vs. \og(5/h) for h = 1/8. 

for each {£, m} are the projection vectors p lm . 

At this point, one may be concerned that the use of monomial basis functions 
might cause G to become ill-conditioned. However, as long as the number of radial 
basis functions is small, this will not be a significant problem. Tests have shown that 
the difference between using the monomial basis and the orthogonal polynomial basis 
for a quadratic fit in r is on the order of the roundoff error expected from solving a 
linear system of that size. 

An additional benefit of this approach comes in cases where the radial derivative 
of $ £m is required at r . This is the case, for example, when the extracted multipole 
amplitudes at ro are to be used as inner boundary data for the evolution of a one- 
dimensional second-order wave equation for each mode as in pQ|. Based on the definition 
of R n in firUj) . it is clear that the coefficient $ 1£m is the derivative (with respect to u) 
of $ £m . Therefore, we have (d$^ m /dr) = & lim /5. Again, this is easily computed 
using projection vectors that are rows of P. 

Fiske's analysis of truncation error in the approximation to the Dirac delta function 
also applies when computing d$ £m /dr, only now it involves odd powers to approximate 
the derivative of the delta function. Here the error will converge as 0(S 2 ) for a linear 
or quadratic fit, and will be 0(<5 4 ) for a cubic fit. Increasing the radial basis to include 
i? 3 will reduce the error in <3> £m by only a very small amount and will not improve the 
convergence properties of that error. However using a cubic fit to compute d$ £m /dr 
offers significant improvement in both the error and its convergence properties, as shown 
in figure |U Here we have a convenient rule of thumb: a quadratic fit in r is more than 
sufficient to compute the multipole amplitudes of $, but a cubic fit is strongly suggested 
with the derivatives of those amplitudes are also required. 
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3-4- Symmetry issues 

One final advantage of the least-squares perspective comes in its application to symmetry- 
issues. Fiske |E] has studied extensively the symmetries of the basis functions Ya and 
shown that the dual basis functions, Y A , have the same symmetries. This is useful 
in cases where the three-dimensional simulation is computed in a single octant with 
explicit reflection symmetries at the coordinate planes. Fiske's analysis shows how to 
apply Misner's method in such cases by exploiting the symmetries of Ya to compute 
the sums used to construct Gab, which must be taken over all octants, using only data 
defined within the computational octant. The result is a block diagonal linear system 
constructed so that coefficients of harmonics that do not share the octant symmetry are 
forced to vanish. 

Again, the least squares perspective offers a slightly different approach to the 
problem. First, because the least squares fit only relies on the basis functions being 
linearly independent, it is possible simply to perform the fit over all spherical harmonics 
using data in the octant. A problem with this approach is that it could cause G to 
become ill-conditioned as the number of basis functions grows large. 

A simple fix is to restrict the least squares fit to only those basis functions that 
share the octant symmetry of the simulation, and ignore all other basis harmonics. For 
example, because the test function in this paper is based on Y 2 q, the function $ will be 
even under reflection across the coordinate plane boundaries of the first octant. Using 
real-valued spherical harmonics, which are proportional to cos(m0) for positive m and 
sin(m0) for negative m, it is easy to show that the only harmonic functions that have 
this octant symmetry are those with even I and even, positive m. This guarantees 
that the fit matrix G and the corresponding linear system are much smaller than those 
which use all spherical basis functions. A fit with £ max = 2, could potentially involve 
nine angular basis functions, but only three of these (loo, ^20, and Y22) share the octant 
symmetry of Y^o- This cuts the size of linear system by a factor of three, significantly 
reducing computation time. While G would no longer as sparse as in Fiske's approach, 
this should not present a significant additional computational burden. 

4. Discussion 

Misner's approach to the multipole decomposition of grid data offers a significant 
computational advantage over the traditional method of first interpolating grid data 
onto the two-sphere, then integrating against the spherical harmonics. For a shell size 
5 that is proportional to the grid spacing h, fourth-order convergence is achieved by a 
simple least-squares fit over a quadratic polynomial in r and the specified set of basis 
functions (spherical harmonics) in {8,<p}. This accuracy is achieved without the need 
to consider coordinate patches that cover the two-sphere. 

The construction of projection vectors that can be computed at the beginning of 
a simulation and stored for use at each computational time level offers a significant 
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increase in computational speed, as well. While such operators can in principle be 
computed for the interpolation/integration method, their construction is significantly 
easier in Misner's method because the projection operators are simply the result of 
solving the linear system associated with a least squares regression problem. 

Re-casting Misner's method as a local least squares problem does offer additional 
advantages in accuracy and computational simplicity. Tests suggest that a (r- 
dependent) tricube weighting scheme gives good results. 

This analysis has not investigated the application of Misner's method to more 
complex problems that would use, for example, tensor harmonics. However, because 
the method is based only on the orthonormality of the basis functions (indeed, on the 
weaker requirement of linear independence) on the two-sphere, it should extend without 
difficulty to those problems. 

We have also not investigated the application of Misner's method to mesh 
refinement. The analysis of Fiske, et al demonstrates the simplicity of Misner's 
method to cases of fixed mesh refinement. Although it is not clear from that paper, a 
strict application of Misner's weighting scheme is slightly more complicated for refined 
grids than for the fixed grids considered here. This is because one must account for the 
fact that points in a refined region belong to cells with smaller volume h 3 . However, 
from the local least-squares perspective the weight of each point only depends on its 
proximity to r$. It is not clear that any advantage is gained by adjusting the weight of 
a point based on the density of grid points nearby. 

Application to adaptive mesh refinement can complicate the process because there 
is no way a priori to establish the distribution of grid points at any given time level. 
This would seem to rule out pre-computing the projection vectors and slow down 
the procedure. However, because the computation is done over a shell that scales 
quadratically with resolution, this may not represent a significant time penalty. 

More importantly, it may not even be necessary to consider adaptive mesh 
refinement when applying Misner's method. As long as the grid size, number of basis 
functions, and grid resolution are sufficient to compute the multipole amplitudes at the 
most coarse level, the multipole decomposition can be applied only to grid points on 
that level, and these are known a priori. Unless high-frequency modes are of interest, 
the values of the wave on the coarse grid should be sufficient to compute the multipole 
amplitudes. 
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